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y—^ Abstract 

^^ , A broad range of inverse problems can be abstracted into the problem of minimizing the 

t"*"- ' sum of several convex functions in a Hilbert space. We propose a proximal decomposition algo- 

^5 , rithm for solving this problem with an arbitrary number of nonsmooth functions and establish 

its weak convergence. The algorithm fully decomposes the problem in that it involves each 
function individually via its own proximity operator. A significant improvement over the meth- 
ods currently in use in the area of inverse problems is that it is not limited to two nonsmooth 
functions. Numerical applications to signal and image processing problems are demonstrated. 
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1 Introduction 

Throughout this paper, 7^ is a real Hilbert space with scalar product (• | •), norm || • ||, and distance 
d. Moreover, {fi)i<i<m are proper lower semicontinuous convex functions from H to ]— oo,+oo]. 
We consider inverse problems that can be formulated as decomposed optimization problems of the 
form 

m 

minimize 2,fi{^)- (1-1) 

In this flexible variational formulation, each potential function /j may represent a prior constraint 
on the ideal solution x or on the data acquisition model. The purpose of this paper is to propose 
a decomposition method that, under rather general conditions, will provide solutions to (1.1). 

To place our investigation in perspective, let us review some important special cases of (1.1) for 
which globally convergent numerical methods are available. These examples encompass a variety 
of inverse problems in areas such as signal denoising [25, 44], signal deconvolution [17], Bayesian 
image recovery [16], intensity-modulated radiation therapy [10, 13], image restoration [5, 6, 15], 
linear inverse problems with sparsity constraints [24, 29, 32, 48], signal reconstruction from Fourier 
phase information [37], and tomographic reconstruction [2, 10, 46]. 

(a) If the functions {fi)i<i<m are the indicator functions (see (2.1)) of closed convex sets 
{Ci)i<i<m in Ti-, (1.1) reduces to the convex feasibility problem [10, 13, 18, 46, 50] 



findxGp|C7i, (1.2) 



i=l 



which can be solved by projection techniques, e.g., [4, 12, 19, 36]. 

(b) The constraint sets in (a) are based on information or measurements that can be inaccurate. 
As a result, the feasibility set fli^i Ci may turn out to be empty. An approximate solution 
can be obtained by setting, for every i G {1, . . . ,m}, fi = oJid^,, where dc^ is the distance 
function to Ci (see (2.2)) and where uji G ]0, 1]. Thus, (1.1) becomes 

m 

minimize 2_,^idci{^)- (1-3) 

This approach is proposed in [17], where it is solved by a parallel projection method. Finite- 
dimensional variants based on Bregman distances are investigated in [11]. 

(c) If the functions {fi)i<i<m-i are the indicator functions of closed convex sets (Ci)i<j<m-i 
in TC and fm- x i-^ W^ ~ ''IP for some r ^ Ti, then (1.1) reduces to the best approximation 



problem [2, 21] 



minimize ||x — r|| . (1-4) 



TTl — 1 



x& f] Ci 



i = l 



Several algorithms are available to solve this problem [7, 21, 33, 34, 49]. There are also 
methods that are applicable in the presence of a more general strictly convex potential fm] 
see [20] and the references therein. 



(d) In [26], the special instance of (1.1) in which m = 2 and /2 is Lipschitz-differentiable on 
7i is shown to cover a variety of seemingly unrelated inverse problem formulations such 
as Fourier regularization problems, constrained least-squares problems, split feasibility prob- 
lems, multiresolution sparse regularization problems, geometry/texture image decomposition 
problems, hard-constrained inconsistent feasibility problems, as well as certain maximum a 
posteriori problems (see also [6, 8, 9, 16, 24, 29, 32] for further developments within this 
framework). The forward-backward splitting algorithm proposed in [26] is governed by the 
updating rule 



Xn+l =Xn + A„f prox^^^j^ (x„ - 7n(V/2(x„) + 6„)) + a„ - Xn) , (1.5) 

where A^ G ]0, 1] and 7„ E ]0, +oo[, where 



1 



P^°^7n/i • ^ ^ argmin 7n/i(y) + ::\\x - y\\ (1.6) 



yen 



2' 



is the proximity operator of 7n/i, and where the vectors a„ and 5„ model tolerances in the 
implementation of prox j-^ and V/2, respectively. Naturally, this 2-function framework can 
be extended to (1-1) under the severe restriction that the functions {fi)2<i<m be Lipschitz- 
differentiable. Indeed, in this case, /2 = X^^2 /« ^^^^ enjoys this property and it can be used 
in lieu of /2 in (1.5). 

(e) The problem considered in [25] corresponds to ?7i = 2 in (1.1). In other words, the smoothness 
assumption on /2 in (d) is relaxed. The algorithm adopted in [25] is based on the Douglas- 
Rachford splitting method [22, 38] and operates via the updating rule 



y„+ 1 = Prox^/2 Vn+an 

Vn+i =yn + K[ prox^j^ (2y„+ 1 - Vn) + bn - y„+ 



where A„ € ]0, 2[ and 7 € ]0,+oo[, and where the vectors a„ and 6„ model tolerances in 
the implementation of the proximity operators. Under suitable assumptions, the sequence 
{yn)n€N converges weakly to a point y & TC and prox^j;^ V ^ Argmin fi + f2- In this approach, 
the smoothness assumption made on /2 in (d) is replaced by the practical assumption that 
prox r^ be implementable (to within some error). 

Some important scenarios are not covered by the above settings, namely the formulations of 
type (1.1) that feature three or more potentials, at least two of which are nonsmooth. In this paper, 
we investigate a reformulation of (1.7) in a product space that allows us to capture instances of 
(1.1) in which none of the functions need be differentiable. The resulting algorithm proceeds by 
decomposition in that each function is involved individually via its own proximity operator. Since 
proximity operators can be implemented for a wide variety of potentials, the proposed framework 
is applicable to a broad array of problems. 

In section 2, we set our notation and provide some background on convex analysis and proxim- 
ity operators. We also obtain closed-form formulas for new examples of proximity operators that 
will be used subsequently. In section 3, we introduce our algorithm and prove its weak convergence. 



Applications to signal and image processing problems are detailed in section 4, where numerical 
results are also provided. These results show that complex nonsmooth variational inverse prob- 
lems, that were beyond the reach of the methods reviewed above, can be decomposed and solved 
efficiently within the proposed framework. Section 5 concludes the paper with some remarks. 



2 Notation and background 

2.1 Convex analysis 

We provide here some basic elements; for proofs and complements see [51] and, for the finite 
dimensional setting, [43]. 

Let C be a nonempty convex subset of Ti. The indicator function of C is 

Jo, ifxGC; 

ic-x^ < (2.1) 

[+00, \i X fC, 

its distance function is 



its support function is 



and its conical hull is 



dc: W^ [0,+oo[: xh^ inf ||x-y||, (2.2) 

y&C 
(Tc: "^^ ^ ] — C)0, +oo] : u I— > sup (x I n), (2.3) 

coneC= J {Ax | x e C}. (2.4) 

A>0 

Moreover, span C denotes the span of C and spanC the closure of span C. The strong relative 
interior of C is 

sri C = {x G C I cone(C — x) = span (C — x) } (2.5) 

and its relative interior is 

ri C = {x G C I cone(C — x) = span (C — x) } . (2-6) 

We have 

intCCsriCCriCcC. (2.7) 

Lemma 2.1 [43, Section 6] Suppose that H is finite- dimensional, and let C and D be convex 
subsets ofTi. Then the following hold. 

(i) Suppose that C ^ 0. Then sriC = riC 7^ 0. 
(ii) vi{C- D) = riC-riD. 
(iii) Suppose that D is an affine subspace and that (riC) DD ^ 0. Then ri(C CiD) = (riC) n D. 



Now let C be a nonempty closed and convex subset of 7i. The projection of a point x m. Ti 
onto C is the unique point Pqx in C such that ||x — Pcx\\ = dc{x). We have 

(Vx G W)(VpG W) p = Pea; 44> [p G C and (Vy G C) {y - p\ x - p) <{)]. (2.8) 

Moreover, dc is Frechet differentiable on TC \ C and 

(VxGHxC) Vcic(x) = ^:P^. (2.9) 

The domain of a function f : 7i ^ ]—oo, +00] is dom/ = {x G 7Y | /(x) < +00} and its set of 
global minimizers is denoted by Argmin /; if / possesses a unique global minimizer, it is denoted by 
argmin g^ f{y)- The class of lower semicontinuous convex functions from TC to ]— 00, +00] which 
are proper (i.e., with nonempty domain) is denoted by ro(7^). Now let / G ro(7^). The conjugate 
of / is the function /* G ro(7^) defined hy f* : TC —<■]— 00, +00] : u 1-^ sup^.^^ {x \ u) — /(x), and 
the subdifferential of / is the set-valued operator 

df:n^2^: x^ {u£n\ (Vy G dom/) {y - x \ u) + f{x) <f{y)]. (2.10) 

We have 

(VxGW) xG Argmin/ <^ G ^/(x) (2.11) 

and 

(^ ^ -wwvy ^ 'w^ /-^(^^ + ■^*(^) - ^^ 1 "^ (1 1o^ 

(VxGW(Vug7^ 1 ./ X , «/ n / , \ Q./ N (2-12) 

yf{x) + f (u) = {x \u) ^ u e df{x). 

Moreover, if / is Gateaux-differentiable at x G 7Y, then df{x) = {V/(x)}. 

Lemma 2.2 Let C be a nonempty closed convex subset of TC, let (p: M ^ M be an even convex 
function, and set f = <j) o dc- Then f G Tq{TC) and f* = gq + </>* o || • ||. 

Proof. Since (/>: M ^ M is convex and even, it is continuous and increasing on [0,+oo[. On the 
other hand, since C is convex, dc is convex. Hence, (p o dc ^s a finite continuous convex function, 
which shows that / G To{TC). Moreover, (j) ode = (Pi^^^yeC II ■ — y||) = infj^GC 4'°\\' ~y\\- Therefore, 

(Vu G TC) f*{u) = sup {x \ u) — inf (j){\\x — y\\) 

= sup {y \ u) + sup {x — y \ u) — {(j) o \\ ■ ||)(x — y) 
yGC xen 

= sup(y I u) + ((/)o II • ||)*(n) 

y&C 

= ac(n) + (</.o||.||)*(u). (2.13) 

Since (</) o || • ||)* = 0* o || • || [31, Proposition 1.4.2], the proof is complete. D 



2.2 Proximity operators 

For detailed accounts of the theory of proximity operators, see [26, Section 2] and [39]. 

The proximity operator of a function / G Tq(7{) is the operator proxj : TC ^ ?{ which maps 
every x ^Ti io the unique minimizer of the function / + ||x — •||^/2, i.e., 

Nx^TL) prox f X = argmin f{y)-\ — ||x — y||^. (2-14) 

yen 2 

We have 

(Vx G 'H){yp G 7Y) p = proxj x <^ x — p & df{p). (2-15) 

In other words, proxj = (Id +df)^^. 

Lemma 2.3 Let f G ro(7^). Then the following hold. 

(i) (Vx G 'H){\/y G TC) \\ proxj x — proxr y||^ < {x — y \ proxr x — proxr y). 
(ii) (Vx G H)(V7 G ]0, +oo[) x = prox^j x + 7proxj,/^(x/7). 

Lemma 2.4 [25, Proposition 11] Let Q be a real Hilbert space, let f G ro(^), and let L: TC ^ Q 
be a bounded linear operator such that L o L* = Kid , for some k G ]0, +oo[. Then / o L G Tq(7{) 
and 

proxjoj^ = Id H — L* o (prox^j — Id ) o L. (2.16) 



2.3 Examples of proximity operators 

Closed-form formulas for various proximity operators are provided in [16, 24, 25, 26, 39]. The 
following examples will be of immediate use subsequently. 

Proposition 2.5 [16, Proposition 2.10 and Remark 3.2(ii)] Set 

f-.n^ ]-oo,+oo] : X ^ ^(/>fc((x I efe)), (2.17) 

fceK 

where: 

(i) 0/KcN; 

(ii) {ckjkeK is an orthonormal basis ofTi; 
(iii) {4>k)k&^ o-re functions in ro(M); 

(iv) Either M. is finite, or there exists a subset L o/ IK such that: 
(a) K \ L is finite; 



(b) (VA; € L) ,/.fe > MO) = 0. 
Then f G TqITI) and 



(Vx E H) proxjx = ^ (prox^^ (x I efc))efc. (2.18) 

fceK 

We shall also require the following results, which appear to be new. 

Proposition 2.6 Let {G, \\ ■ ||) be a real Hilhert space, let L: 7i ^> Q he linear and hounded, let 
z £ G, let J G ]0, +oo[, and set f = j\\L ■ — z|p/2. Then f G Tq{TC) and 

(VxeW) pioxj-x = (ld+-fL*L)-^{x + jL*z). (2.19) 

Proof. It is clear that / is a finite continuous convex function. Now, take x and p in TC. Then 
(2.15) yields p = proxj x 4^ x — p = V(7||L • — z|p/2)(p) <^ x — p = ^L*{Lp — z) 4^ p = 
(Id+7L*L)-i(x + 7L*z). D 

Proposition 2.7 Let C he a nonempty closed convex suhset ofTi, let (p: M — s- M he an even convex 
function which is differentiahle onM.\ {0}, and set f = (f)o dc- Then 



(Vx S Ti) prox /■ X = < 



pmx^,dc[x) -f ^ ( \^ n^^m 

X H — — [Pcx - x), if dc{x) > max 50(0); 

dcia;) (2.20) 

^P^x, if dc{x) < max 90(0). 



Proof. As seen in Lemma 2.2, / € ro('H). Now let x G W and set p = prox^-x. Since is a finite 
even convex function, 90(0) = [— /3, /3] for some /? G [0,+oo[ [43, Theorem 23.4]. We consider two 
alternatives. 

(a) p E C: Let y G C. Then f{y) = 4>[dc{y)) = 0(0) and, in particular, f{p) = 0(0). Hence, it 
follows from (2.15) and (2.10) that 

(y-p|x-p) + 0(O) = (y-p|x-p) + /(p)</(y) = 0(O). (2.21) 

Consequently, {y — p\x— p)<Q and, in view of (2.8), we get p = Pcx. Thus, 

p G C -^ p = Pcx. (2.22) 

Now, let u G df{p). Since p G C, (ic(p) = and, by (2.3), ac{u) > {p\u). Hence, (2.12) 
and Lemma 2.2 yield 

- ||u|| = < aciu) -{p\u)= ac{u) - f{p) - f*{u) = -0(0) - 0*(|ln||). (2.23) 

We therefore deduce from (2.12) that ||n|| G 90(0). Thus, u G df{p) =^ \\u\\ < (3. Since (2.15) 
asserts that x—p G df{p), we obtain ||x— p|| < (3 and hence, since p G C, dc{x) < ||x— p|| < (3. 
As a result, 

p£C => dc{x) < 13. (2.24) 



(b) p ^ C: Since C is closed, dc{p) > and cj) is therefore differentiable at dc{p)- It follows from 
(2.15), the Frechet chain rule, and (2.9) that 

x-p = np) = ^^^ip-Pcp). (2.25) 

Since (/>' > on ]0, +oo[, upon taking the norm, we obtain 

\\p-x\\ = ^'{dc(j})) (2.26) 

and therefore 

P-X= "^~f (PcP-p). (2.27) 

dc{p) 

In turn, appealing to Lemma 2.3(i) (with / = lc) and (2.8), we obtain 

\\Pcp - Pcxf <{p-x\Pcp- Pcx) = ^^^^{PcP -p\PcP- Pcx) < 0, (2.28) 

dc{P) 

from which we deduce that 

Pcp = Pcx. (2.29) 

Hence, (2.27) becomes 



P-^= II ^ j^\, {Pcx-p), (2.30) 

\\P- Pcx\\ 



which can be rewritten as 



Taking the norm yields 



p — X = T, r, r, ;; — w(Pcx — x). (2.31) 

\\p - x\\ + \\p - Pcxr "^ ' "- ' 



\\P - ^11 = T\ if^ii''" p n dcix), (2.32) 

\\p- x\\ + \\p- Pcx\\ 



and it follows from (2.29) that 

dc{x) = \\p - x\\ + \\p - Pcx\\ = \\p - x\\ + dc{p)- (2.33) 

Therefore, in the light of (2.26), we obtain 

dc{x) - dc{p) = \\p - x\\ = 4>'{dc{p)) (2.34) 

and we derive from (2.15) that 

dc{p) = prox^ (ic(x). (2.35) 

Thus, Lemma 2.3(ii) yields 

dc{x) - dc{p) = dc{x) - prox^ dc{x) = prox^, dc{x) (2.36) 

and, in turn, (2.34) results in 

\\p - x\\ = dc{x) - dc{p) = prox^, dc{x). (2.37) 



To sum up, coming back to (2.31) and invoking (2.33) and (2.37), we obtain 

Wd — x\\ 

p^C =^ P = x + - ' , ii' II {Pcx - x) 

\\p- x\\ + \\p- Pcx\\ 

prox^.dc(x) 

= ^^ lh—\ [Pcx-x). (2.38) 

dc{x) 

Furthermore, we derive from (2.35) and (2.15) that 

piC ^ dc{p) > ^ prox^(ic(x) / ^ dc{x) i 90(0) ^ dc{x) > (3. (2.39) 

Upon combining (2.24) and (2.39), we obtain 

peC ^ dc{x) < (5. (2.40) 

Altogether, (2.20) follows from (2.22), (2.38), and (2.40). D 

The above proposition shows that a nice feature of the proximity operator oi (p o dc is that it 
can be decomposed in terms of prox ., and Pq- Here is an apphcation of this result. 



Proposition 2.8 Let C he a nonempty closed convex subset ofTi, let a € ]0, +oo[, letp G [1, +oo[, 



and set f = ad^. Then the following hold. 



(i) Suppose that p = 1. Then 



{a 
x + -—— -{Pcx-x), if dcix) > a; 
dcix) (2.41) 

Pcx, if dc{x) < a. 

(ii) Suppose that p > 1. Then 

\x+ ^^,'A Pcx-x), ifx4C; , ^ 

(VxGW) prox^x=<^ dc{xy "" 7' ^ ^ ' ^2.42) 

yx, if X € C, 

where u[x) is the unique real number in [0, +cx)[ that satisfies i^(x) + {i'[x)/{ap))^'^^~^' = 
dcix). 

Proof (i): Set cj) = a\ ■ \. Then max 00(0) = max [—a, a] = a and (p* = i[-a,a]- Therefore, 
prox^* = P[-a,a] and hence (Vfi G ]a,+oo[) prox^, fi = a. In view of (2.20), we obtain (2.41). 

(ii): Let x € TC and note that, since C is closed, dc{x) > ^ x ^ C. Now set cj) = a\- 1^. Then 
max (90(0) = max{0} = and 0* : ^ i-^ {p-l){ap)^^^^~P^\id\P/^P~^^ /p. Hence, it follows from (2.15) 
and [24, Corollary 2.5] that prox^* dcix) is the unique solution i>{x) € [0, +oo[ to the equation 
dc{x) - v{x) = (t)*'{iy{x)) = (z/(x)/(ap))i/(P-i). Appealing to (2.20), we obtain (2.42). D 

Let us note that explicit expressions can be obtained for several values of p in Proposition 2.8(ii). 
Here is an example that will be used subsequently. 

9 



Example 2.9 Let C be a nonempty closed convex subset of 7i, let a G ]0, +cx3[, and set / = adj . 



Then 



(Vx G Ti) proxj x = < 



, 9a2(Vl + 16dc(x)/(9a2)-l) 
" + M^) iPcx-xl .ixiC- ^2.43) 



X, 



if X G C. 



Proof. Set p = 3/2 in Proposition 2.8(ii). D 



3 Algorithm and convergence 

The main algorithm is presented in section 3.1. In section 3.2, we revisit the Douglas- Rachford 
algorithm in the context of minimization problems (Proposition 3.2), with special emphasis on its 
convergence in a specific case (Proposition 3.3). These results are transcribed in a product space 
in section 3.3 to prove the weak convergence of Algorithm 3.1. 

3.1 Algorithm 

We propose the following proximal method to solve (1.1). In this splitting algorithm, each function 
fi is used separately by means of its own proximity operator. 



Algorithm 3.1 For every i G {1, . . . ,?ti}, let {ai^n)nen be a sequence in Ti. A sequence {xn)n&i 
is generated by the following routine. 

Initialization 

7 G]0,+cx)[ 

(wi)i<i<m G ]0, 1]"^ satisfy YZi ^i = 1 

(yi,o)l<J<m G "^"^ 



Xq 



i=l 

For n = 0, 1, . . . 



,m 



For i = 1, 

[ Pi^n = prox^y^/^^ yi^n + fli.r 



Pn = ^ (^iPi,n 



i=l 

A„G]0,2[ 

For i = 1, . . . ,m 

|_ yj,n+l — yi,n ~r ■^n\^Pn Xn Pi,n) 
Xn+1 — Xfi -\- '^n\Pn Xn)- 



(3.1) 



10 



At iteration n, the proximal vectors {pi,n)i<i<mi as well as the auxiliary vectors {yi,n)i<i<m, 
can be computed simultaneously, hence the parallel structure of Algorithm 3.1. Another feature of 
the algorithm is that some error aj^„ is tolerated in the computation of the ith proximity operator. 



3.2 The Douglas-Rachford algorithm for minimization problems 

To ease our presentation, we introduce in this section a second real Hilbert space {Ti, ||| • |||). As 
usual, ^ denotes weak convergence. 

The (nonlinear) Douglas-Rachford splitting method was initially developed for the problem of 
finding a zero of the sum of two maximal monotone operators in [38] (see [22] for recent refine- 
ments). In the case when the maximal monotone operators are subdifferentials, it provides an 
algorithm for minimizing the sum of two convex functions. In this section, we develop this point 
of view, starting with the following result. 

Proposition 3.2 Let fi and /2 be functions in TQ{7i), let (a„)„gN (^nd {bn)neN be sequences in 
Ti, and let (?/„)ngN be a sequence generated by the following routine. 

Initialization 
7 G]0,-Foo[ 
2/0 eTi 

For n = 0, 1, . . . (3.2) 

A„, G]0,2[ 

Vn+l =yn + KyP^OX^f^ ^'^yn+\ " Vn) + K - 2/„+i 

Set 

G = ATgnunfi+ f2 and T = 2 prox^^.^ o (2 prox^_^_^ - Id ) - 2 prox^^^ -|- Id , (3.3) 

and suppose that the following hold. 

W ,„ m™ fi(.^) + f2{x) = +oo. 
(ii) € sri(dom/]^ — dom/2). 
(iii) EneN -^n(2 - An) = +00. 

(iv) E„GNAn(|||an||| + |||&n|||) < +00. 

Then G ^ 0, (y„)neN converges weakly to a fixed point y of T, and prox^f y G G. 



11 



Proof. It follows from (ii) that doui{fi + /2) = douifi D doin/2 7^ 0. Hence, since fi + /2 is 
lower semicontinuous and convex as the sum of two such functions, we have /i + /2 ^ ro(7i). In 
turn, we derive from (i) and [51, Theorem 2.5.1(ii)] that 



G/0. 



(3.4) 



Next, let us set Ai = dfi, A2 = 9/2, and Z = {xeTC\Oe Aix + A2X}. Then Ai and A2 are 
maximal monotone operators [51, Theorem 3.1.11]. In addition, in view of (2.15), the resolvents 
of 7A1 and 7A2 are respectively 



J-yAi = (Id +7A1) = prox^^^ and J^Aj = (Id +7A2) 
Thus, the iteration in (3.2) can be rewritten as 

A„G]0,2[ 



prox^^^ 



(3.5) 



(3.6) 



Vn) +K- ?/„+ 

Moreover, it follows from (2.11), (ii), and [51, Theorem 2.8.3] that 

G = {a; G ?i I G d{f^ + f2){x)] = {a; G ?i | G df^{x) + df^ix)] = Z. 



(3.7) 



Thus, (3.4) yields Z ^ and it follows from (iii), (iv), and the results of [22, Section 5] that 
(yn)neN converges weakly to a fixed point y of the operator 2J^Ai ° (2 J7A2 — Id ) — 2J^A2 + Id > 
and that J7A2I/ ^ ■^- In view of (3.3), (3.5), and (3.7), the proof is complete. D 

It is important to stress that algorithm (3.2) provides a minimizer indirectly: the sequence 
(yn)neN is first Constructed, and then a minimizer of fi + /2 is obtained as the image of the weak 
limit y of {yn)neN under prox j. . In general, nothing is known about the weak convergence of 
the sequences (prox^^ Z/n)neN and (prox^^, y„)neN- The following result describes a remarkable 
situation in which (prox j^ Z/n)neN does converges weakly and its weak limit turns out to be a 
minimizer of fi + f2- 

Proposition 3.3 Let D be a closed vector subspace ofTi, let f G Tq(T-C), let (a„)„gp^ be a sequence 
in Ti, and let {xn)nefi ^6 a sequence generated by the following routine. 

Initialization 

7 G]0,+oo[ 

xo = Pd yo 

For n = 0, 1, . . . 



y„+ 1 = prox^^ y„ + a„ 
Pn = Pd yn+^ 
A„G]0,2[ 

Vn+l — yn~^ ^n (2p„, — Xn " 
Xn-\-l — Xn + '^n\Pn Xn). 



(3.8) 



Vr 



Let G be the set of minimizers of f over D and suppose that the following hold. 
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(i) lim f{x) = +00. 

xSiD, |||a;|||— »+oo 

(ii) G sn{D-domf). 

(iii) EnGN -^n(2 " A„) = +00. 
(iv) EnGN-^nlll«nIII < +°°- 

Then G 7^ and {xn)neN converges weakly to a point in G. 



Proof. Set /i = LD, /2 = /, and (Vn G N) b„ = 0. Then (2.1) and (2.14) yield prox^/^ = Pd 
and, since D is a closed vector subspace, Pd is a linear operator. Hence, proceeding by induction, 
we can rewrite the update equation for a;„ in (3.8) as 

= PDyn + K {^Pd Pn - PoXn- Pd ^n+i) 

= Pd (Vn + K{2Pn -Xn- yn+\) ^ 
= PDyn+l- 

As a result, (3.8) is equivalent to 

Initialization 

7 G]0,+oo[ 

For n = 0, 1, . . . 

Xn = Pd yn 

y„+i =prox^^y„ + a„ 



(3.9) 



(3.10) 



Pn 



PdVt 



AnG]0,2[ 

yn+1 =yn + >'n{2Pn 



Xr 



Vr 



Thus, since 



(yxen)iyyen) PD{2y-x)=2PDy-PDX, (3.11) 

(3.10) appears as a special case of (3.2) in which we have introduced the auxiliary variables a;„ 
and p^. In addition, the operator T of (3.3) becomes 

T = 4(Pd o prox^^f ) -2Pd-2 prox^^^ + Id . (3.12) 

Since (i)-(iv) are specializations of their respective counterparts in Proposition 3.2, it follows from 
Proposition 3.2 that G ^ and that there exists a fixed point y of T such that yn ^ y and 
prox^^ y G G. Note that, since G C D, prox^^y G D and, in turn, PD(prox^^y) = prox^^j/. 
Thus, in view of (3.12), we obtain 



Ty = y <^ 4Pr)(prox^^ y) - 2PDy - 2 prox^^ y + y 
•^ 2PD(prox^^ y) - Poy = prox^^ y 
•^ prox^f y = Poy- 



y 



(3.13) 
(3.14) 
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Hence, since prox ^ y (^ G, we also have PdV G G- On the other hand, since Pd is hnear and 
continuous, it is weakly continuous and therefore y^ ^ y ^ PdVu ~^ PdV £ G. We conclude 
that Xn -^ PdV € G. D 



3.3 Convergence of Algorithm 3.1 

The convergence of the main algorithm can now be demonstrated. 

Theorem 3.4 Let G he the set of solutions to (1.1) and let {xn)n&i be a sequence generated by 
Algorithm 3.1 under the following assumptions. 

(i) lim fl{x)^ V fm.{x) = +00. 

j|a;||— >+oo 

(ii) (0, . . . , 0) € sri { (x - xi, . . . , X - Xm) | x G ?^, xi e dom /i, . . . , x^ G dom /„} . 

(iii) EnGN -^n(2 " A„) = +00. 

(iv) (Vi e {l,...,m}) EnGN-^rillai,™!! < +00. 

Then G ^ and (x„)„i=n converges weakly to a point in G. 

Proof. Let 1-L be the real Hilbert space obtained by endowing the TTi-fold Cartesian product Ti"^ 
with the scalar product 

m 

{{■\-))- {x,y)^^i^i{xi\yi), (3.15) 

i=l 

where (a;j)i<i<m is defined in (3.1), and where x = (xj)i<j<m and y = {yi)i<i<m denote generic 
elements in Ti. The associated norm is denoted by ||| • |||, i.e.. 



: X 



^(Ji||xj||2. (3.16) 



Furthermore, set 
and 



D = {{x,...,x) e7i\ xen] (3.17) 



/:?i^]-oo,+oo]:a;^^/i(xi). (3.18) 



It follows from (3.16) that D is a closed vector subspace of Ti with projector 



^ m m \ 

Pd: x^ f ^ ujiXi, • • • , ^ uJiXi j . 



(3.19) 



and that the operator 

j-.n^D: x^{x,...,x) (3.20) 
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is an isomorphism. In addition, / € ro(7i) and we derive from (2.14), (3.16), and (3.18) that 

prox_^ : x^ { prox^j/^^ xi, . . . , prox^^/^^ Xm) ■ (3.21) 

From the sequences {xn)nen, (Pn)neN, {{yi,n)nen)i<i<m, {{Pi,n)nm)i<i<m, and ((aj,„)„gN)i<i<m of 
Algorithm 3.1 we define, for every n € N, 



X 



n =j{Xn), Pn = J{Pn), Vn = {yi,n)l<i<m, 2/„+l/2 = {Pi,n)l<i<m, and a„ = {ai^n)l<i<m- (3.22) 



It follows from (3.19), (3.20), and (3.21) that the sequences defined in (3.22) are precisely those 
involved in (3.8), and that the set of minimizers G in Proposition 3.3 is precisely 

G = j{G). (3.23) 

On the other hand, it follows from (3.16), (3.17), and (3.18) that the properties (i)-(iv) above yield 
their respective counterparts in Proposition 3.3. Thus, we deduce from Proposition 3.3 and (3.23) 
that {xn)neN converges weakly to a point j{x) for some x € G. Thus, {xn)neN = {j^^{xn))neN 
converges weakly to x € G. D 

Remark 3.5 

(i) We have conveniently obtained Algorithm 3.1 as a direct transcription of a special case (see 
Proposition 3.3) of the Douglas- Rachford algorithm transposed in a product space. A similar 
decomposition method could be obtained by using the theory of partial inverses for monotone 
operators [45]. 

(ii) When m = 2, Algorithm 3.1 does not revert to the standard Douglas- Rachford iteration 
(1.7). Actually, even in this case, it seems better to use the former to the extent that, as 
seen in Theorem 3.4, it produces directly a sequence that converges weakly to a minimizer 

0f/l + /2. 

To conclude this section, we describe some situations in which condition (ii) in Theorem 3.4 is 
satisfied. 

Proposition 3.6 Set C = {{x — xi, . . . ,x — x-m) \ x G "H, xi S dom/i, . . . , x^ G dom/m} and 
suppose that any of the following holds. 

(i) C is a closed vector subspace. 

(ii) fXiLidomfi ^ and (dom/j)i<j<m are affine subspaces of finite dimensions. 

(iii) f]^^douifi 7^ and (dom/j)i<j<m are closed affine subspaces of finite codimensions. 

(iv) G intC. 

(v) dom /i n ni^2 i^t do™ /j 7^ 0- 

(vi) 7i is finite- dimensional and fX^LiTidomfi ^ 0. 
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Then € sriC. 

Proof. We use the notation of the proof of Theorem 3.4, hence C = D — doni/. 

(i): We have C = span C. Since C C coneC C span C C spanC, we therefore obtain 
coneC = span C. Appeahng to (2.5), we conclude that € sriC. 

(ii)^(i): The assumption imphes that dom/ = dom/i x • • • x douifm is a finite-dimensional 
affine subspace of 7i and that D n dom/ ^ 0. Since D is closed vector subspace, it follows from 
[30, Lemma 9.36] that D — dom/ is a closed vector subspace. 

(iii)=^(i): Here dom/ = dom/i x ••• x dom/m is a closed affine subspace of 7i of finite 
codimension and that D f] dom/ ^ 0. Appealing to [30, Theorem 9.35 and Corollary 9.37], we 
conclude that D — dom / is a closed vector subspace. 

(iv): See (2.7). 

(v)^(iv): See the proof of [3, Theorem 6.3]. 

(vi): Using Lemma 2.1(i)&:(ii), we obtain € sriC <^ € sri(D — dom/) <^ € ri(Z) — dom/) 
44> € ri Z) — ri dom f = D — t'\ dom f <^ D Dri dom / ^ <^ Cth=i ^i dom fi ^ 0. D 



4 Applications to signal and image processing 

To illustrate the versatility of the proposed framework, we present three applications in signal and 
image processing. In each experiment. Algorithm 3.1 is implemented with cvi = 1/m, Xn = 1.5, and, 
since the proximity operators required by the algorithm will be computable in closed form, we can 
dispense with errors and set aj^„ = in (3.1). As a result, conditions (iii) and (iv) in Theorem 3.4 
are straightforwardly satisfied. In each experiment, the number of iterations of the algorithm is 
chosen large enough so that no significant improvement is gained by letting the algorithm run 
further. 



4.1 Experiment 1 

This first experiment is an image restoration problem in the standard Euclidean space Ti. = 
where A'^ = 512. The original vignetted N x N image x is shown in figure 1 (the vignetting is 
modeled by a black area in the image corners) . The degraded image z shown in figure 2 is obtained 
via the degradation model 

z = Lx + w, (4.1) 

where L is the two-dimensional convolution operator induced by a 15 x 15 uniform kernel, and 
where tt; is a realization of a zero-mean white Gaussian noise. The blurred image-to-noise ratio is 
201og]^Q(||Lx||/||t(;||) = 31.75 dB and the relative quadratic error with respect to the original image 
is 201ogio(||z-x||/||x||) = -19.98 dB. 
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Figure 1: Experiment 1. Original image. 




Figure 2: Experiment 1. Degraded image. 
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Figure 3: Experiment 1. Image restored with 300 iterations of Algorithm 3.1 (7 = 1/4). 



The pixel values are known to fall in the interval [0, 255]. In addition, the vignetting area S of 
the original image is known. This information leads to the constraint set 



Ci = [0,255]^ njx G W I a;l§ = 0}, 



(4.2) 



where rrls denotes the coordinatewise multiplication of x with the characteristic vector 1§ of S (its 
/cth coordinate is 1 or according as /c G § or A; ^ S), and where the zero image. The mean value 
// G ]0, 255 [ of X is also known, which corresponds to the constraint set 



C2 = {xG W I (x I 1) =iVV}, 



(4.3) 



where \ = [1, . . . , 1] G M . In addition, the phase of the discrete Fourier transform of the original 
image is measured over some frequency range B C {0, . . . , A^^ — 1} [17, 37, 42]. If we denote by 
^ = (IXfel exp(iZxfc))Q<^<jy2_;^ the discrete Fourier transform of an image x £ H and by (0fc)fcGiD) 
the known phase values, we obtain the constraint set 



C3 = {xen\iykG D) zxk 

A constrained least-squares formulation of the problem is 

\Lx — z\\ 



J- 



or, equivalently. 



minimize 

xGCinC2nC3 



minimize in-Ax) + \\Lx 
xeCinCi 



(4.4) 

(4.5) 
(4.6) 
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However, in most instances, the phase cannot be measured exactly. This is simulated by intro- 
ducing a 5 % perturbation on each of the phase components {(j)k)k&B- To take these uncertainties 
into account in (4.6), we replace the "hard" potential tc^ by a smoothed version, namely ad^ , for 
some a G ]0, +00 [ and p G [1, +00 [. This leads to the variational problem 

minimize acR, (x) + \\Lx — z\\ , (4.7) 

x€CinC2 ^ 



-z 



2 



which is a special case of (1.1), with ?Ti = 4, /i = lcx-, /2 = i-c-i-, fs = Oid^ , and f^ = \\L 

Let us note that, since Ci is bounded, condition (i) in Theorem 3.4 is satisfied. In addition, 

it follows from Proposition 3.6(vi) that condition (ii) in Theorem 3.4 also holds. Indeed, set 

E = ]0, 255[^' n yl n C2, where A = {x e H \ xlg = O}. Then it follows from (4.3) that 

^'^ a-ls)Gi?. (4.8) 



iV2 _ card S 
Hence, since A and C2 are affine subspaces, (4.2) and Lemma 2.1(iii) yield 

4 

2, 



Pi ri dom fi = ri Ci n ri ^2 = (ri Ci) n ^2 = (ri [0, 255]^ )nAnC2 = Ey^0. (4.9) 



Problem (4.7) is solved for the following scenario: B corresponds to a low frequency band 
including about 80 % of the frequency components, p = 3/2, and a = 10. The proximity operators 
required by Algorithm 3.1 are obtained as follows. First, prox^ and proxj^ are respectively the 
projectors onto Ci and C2, which can be obtained explicitly [18]. Next, proxj- is given in Exam- 
ple 2.9. It involves Pc^, which can be found in [18]. Finally, proxj^ is supplied by Proposition 2.6. 
Note that, since L is a two-dimensional convolutional blur, it can be approximated by a block cir- 
culant matrix and hence (2.19) can be efficiently implemented in the frequency domain via the fast 
Fourier transform [1]. The restored image, shown in figure 3, is much sharper than the degraded 
image z and it achieves a relative quadratic error of —23.25 dB with respect to the original image 

X. 



4.2 Experiment 2 

In image recovery, variational formulations involving total variation [15, 44, 47] or sparsity pro- 
moting potentials [5, 8, 14, 28] are popular. The objective of the present experiment is to show 
that it is possible to employ more sophisticated, hybrid potentials. 

In order to simplify our presentation, we place ourselves in the Hilbert space Q of periodic 
discrete images y = {'r]k,i){k,i)&'i? with horizontal and vertical periods equal to N {N = 512), 
endowed with the standard Euclidean norm 



\ 



N-lN-l 



E E i^mP- (4.10) 



A:=0 «=0 
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As usual, images of size N x N are viewed as elements of this space through periodization [1]. The 
original 8-bit satellite image y (z Q displayed in figure 4 is degraded through the linear model 

z = Ly + w, (4.11) 

where L is the two-dimensional periodic convolution operator with a 7 x 7 uniform kernel, and w 
is a realization of a periodic zero-mean white Gaussian noise. The resulting degraded image z a G 
is shown in figure 5. The blurred image-to-noise ratio is 201ogxo(||-^y||/||u^||) = 20.71 dB and the 
relative quadratic error with respect to the original image is 201ogxo(||-2 — y||/||y||) = —12.02 dB. 

In the spirit of a number of recent investigations (see [16] and the references therein), we use 
a tight frame representation of the images under consideration. This representation is defined 
through a synthesis operator F*, which is a linear operator from TC = M to Q (with K > N'^) 
such that 

F*oF = Kld (4.12) 

for some k G ]0, -|-oo[. Thus, the original image can be written as y = F*x, where x £ Ti is a vector 
of frame coefficients to be estimated. The rationale behind this approach is that, by appropriately 
choosing the frame, a sparse representation x of y can be achieved. 

The restoration problem is posed in the frame coefficient space TY. We use the constraint set 
imposing the range of the pixel values of the original image y, namely 

C = {xen\ F*x£D}, where D = {y e Q \ (V(fc,/) G {0, . . . , N - 1}^) i]k,i G [0,255]}, (4.13) 

as well as three potentials. The first potential is the standard least-squares data fidelity term 
X I— > \\LF*x — z|p. The second potential is the £^ norm, which aims at promoting a sparse 
frame representation [16, 28, 48]. Finally, the third potential is the discrete total variation tv, 
which aims at preserving piecewise smooth areas and sharp edges [15, 44, 47]. Using the notation 
(%,0a,. /-jg^a = (^/,fc){A:,/)ez2; the discrete total variation of y G ^ is defined as 

N-lN-l 

My) = E E ^M(Viy,(Vi(yT))T), (4.14) 

fc=0 1=0 

where Vi : ^ ^ M^^^ is a discrete vertical gradient operator and where, for every {k,l,q,r} C 
{0, ...,A^-1}, we set 

Qk,l = Qk,l,k,l , (4.15) 

with 

ek,l,,,r: M^x^ X M^x^ -. R: (K,4<a,.<Ar-i' [^a,t>l^a,t<N-i) ^ ^/M + M ■ (4-16) 

A common choice for the gradient operator is Vi : y i-^ [%+!,« — Vk,i]o<k,i<N~i- As is customary in 
image processing [35, Section 9.4], we adopt here a horizontally smoothed version of this operator, 
namely, 



1 



Vi: ^^M^xA^. y^ ^^^^^^^^_^^^^^^^^^^^_rj^.^^-j^^^^^^_^_ (4.17) 
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Figure 4: Experiment 2. Original image. 



We thus arrive at a variational formulation of the form (1.1), namely 

minimize ic{^) ~^ \\LF*x — z\\ + a|[a;||^i + /?tv(F*x), 



(4.18) 



where a and f3 are in ]0, +co[. Since C is bounded, condition (i) in Theorem 3.4 is satisfied. In 
addition, it is clear from Proposition 3.6(vi) that condition (ii) in Theorem 3.4 also holds. Indeed, 
all the potentials in (4.18) have full domain, except lq- However, Lemma 2.1(i) implies that 
ri dom tc* = ri C 7^ since € C. 

Although (4.18) assumes the form of (1.1), it is not directly exploitable by Algorithm 3.1 
because the proximity operator of tvoF* cannot be computed explicitly. To circumvent this 
numerical hurdle, the total variation potential (4.14) is split in four terms and (4.18) is rewritten 
as 



minimize \\LF*x — z\\ +a||x|Li+/3> tv,-(F*x), 

II II II lit f ^_^ v /, 



x&C 



(4.19) 



i=0 



where 



Af/2~1 N/2~l 
■y^ ^ Yl ^2fc+q,2i+r(Viy, (Vi(y^))'^). 



(V(g,r)e{0,l}2) tv,+2r:a 

fe=0 /=0 

For every q and r in {0, 1}, let [g^r be the decimation operator given by 

Tn>2Afx2Ar 



I 



Vqy 
oNxN . 



q,r- 



: V 



[vk. 



0<k,l<2N-l 



[V2k 



+g,2/+rJo<fc;<Af-l' 



(4.20) 



(4.21) 
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Figure 5: Experiment 2. Degraded image. 




Figure 6: Experiment 2. Image restored by (4.27), using 350 iterations of Algorithm 3.1 with 
7 = 150. 



22 




Figure 7: Experiment 2. Image restored without the total variation potential in (4.27), using 350 
iterations of Algorithm 3.1 with 7 = 150. 






Figure 8: Experiment 2. Image restored without the i^ potential in (4.27), using 350 iterations of 
Algorithm 3.1 with 7 = 150. 
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and set 

Uq+2r '■ Q 

where Vi is defined in (4.17), 



and 



Moreover, set 






h: 



oNxN 



oiVxTV 



pTVxAT 



r,NxN 



1 



■ y^ i 



q,r 



Voy Viy 

(Vi(y^))T Vay 



Af/2-1 N/2-1 

{v,v). 



y 



N/2-1 N/2-1 
^ '^ 2^ 2^ Qk,l+N/2,k+N/2,l[ 



k=0 1=0 



Then it follows from (4.20) and (4.22) that 

(Vi e {0,1,2,3}) tv, = hoUi. 
Hence, (4.19) becomes 

3 

minimize \\LF*x — z\\ +a||x|Li+/3> hiUiF* 
„ II II II III. ^_^ \ 



xGC 



(4.22) 

(4.23) 
(4.24) 

(4.25) 
(4.26) 
(4.27) 



j=0 



Problem (4.27) is a specialization of (1.1), in which m = 7, fi = ic^ /2 = H-^-^* •— z|p, /s = a\\ ■ ||^i, 
and /j+4 = (3 hoUioF* for i G {0, 1, 2, 3}. To implement Algorithm 3.1, we need the expressions of 
the proximity operators of these functions. The proximity operator of /i can be calculated by first 
observing that the projection onto the set D of (4.13) is explicit, and by then applying Lemma 2.4, 
which states that (4.12) and (4.13) imply that 

proxj^ = prox,^„^* = Id +-F o (prox,^ - Id ) o F* = Id +-F o {Pj^ - Id) o F* . (4.28) 

On the other hand, the proximity operator of /2 can be derived from Proposition 2.6 using a 
frequency domain implementation (as in section 4.1), and by again invoking Lemma 2.4. Next, the 
proximity operator of /s can be found in [26, Example 2.20]. Finally, the operators (proxj.)4<j<7 
are provided by the following fact. 



Proposition 4.1 SetH: 



oNxN 



xN 



■ ^ = Mo<kA<N-l ^ [^kAoKkAKN-V "^^^^^ 



(V(A:,/)G{0,...,iV/2-l}2) 



'^k,i = J^k,i 

T^k+N/2,l+N/2 = ^k+N/2,l+N/2 
'^k,l+N/2 = (^k.liv) ^k,l+N/2 
^'^k+N/2,1 = (^k,l{v) l'k+N/2,1 



with 



o-k,l ■■ v^ < 



1- 
10, 



K/3 



^k,l+N/2\'^ + Wk+N/2,l\'^ 



, if ^JWk,l+N/2\^ + Wk+N/2,l\^ ^ '^/'; 



(4.29) 



otherwise. 
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Then, for every i € {0, 1, 2, 3}, 

prox^ = Id +-F o{U* onoUi-ld)o F*. (4.30) 

Proof. Set v?: M^ ^ M: (6,^2) ^ K/3y^|^i|2 + |^2p- By applying Proposition 2.8(i) in M^ ^ith the 
set {(0,0)}, we obtain 



(V(ei,e2)GM2) prox^(6,6) = < 



('-7rap)'^"«''' «^^^^^^^^'"'^ (4.3.) 

0, otherwise. 



Now set p = [iTk,i]o<k,i<N-i = prox^^^ti. In view of (2.14), (4.25), and (4.16), p minimizes over 
p G M^x^ the cost 

N/2-1 N/2-1 ^ N-lN-l 



kPHp) + -\\v - pf = kP J2 ^ Qk,l+N/2,k+N/2AP^P) + ^J2J2 



r, ^^ ^_^ \^k,l — T^k,l\ 
k=0 1=0 k=0 1=0 



N/2-1 N/2-1 



Yl Y [l^P^J\T^k,l+N/2\^ + \T^k+N/2,l\'^ 
k=0 1=0 

+ 2 (l^*:,'+^/2 ~ ■'^fc,«+Af/2i + Wk+N/2,l — T^k+N/2,l\ )) 

N/2-1 N/2-1 

"*" 2 ^ ^ (kA:,/-*fe,d +Wk+N/2,l+N/2-T^k+N/2,l+N/2\ )■ (4.32) 
fc=0 1=0 



Therefore, 



(V(fc,0G{0,...,iV/2-l}2) <^ 



{T^k,l+N/2j'^k+N/2,l) = PI'OX^(l/fe^i+7v/2; 2^fc+Af/2,0) 

vTfc.i = i^fc,i, (4.33) 

^^k+N/2,l+N/2 = '^k+N/2,l+N/2- 



Appealing to (4.29) and (4.31), we obtain 11 = prox^^/j. Now, let i € {0,1,2,3}. It follows 
from (4.22) that Ui is a separable two-dimensional Haar-like orthogonal operator [35, Section 5.9]. 
Hence, appealing to (4.12), we obtain {Ui o F*) o (U^ o F*)* = kH . In turn. Lemma 2.4 yields 

prox/^+4 =Prox^fco{c/,oF*) 

= Id +-{Ui o F*y o (prox,^^ - Id ) o ([/, o F*) 

= ld+-Fo(U*oUo Ui-ld)oF*, (4.34) 

K 

which completes the proof. D 

In (4.27), we employ a tight frame (k = 4) resulting from the concatenation of four shifted 
separable dyadic orthonormal wavelet decompositions [41] carried out over 4 resolution levels. The 
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shift parameters are (0,0), (1,0), (0,1), and (1,1). In addition, symlet filters [27] of length 8 are 
used. The parameters a and /? have been adjusted so as to minimize the error with respect to 
the original image y. The restored image we obtain is displayed in figure 6. It achieves a relative 
mean-square error with respect to y of —14.82 dB. For comparison, the result obtained without the 
total variation potential in (4.27) is shown in figure 7 (error of —14.06 dB), and the result obtained 
without the £^ potential in (4.27) is shown in figure 8 (error of —13.70 dB). It can be observed 
that the image of figure 7 suffers from small visual artifacts, whereas the details in figure 8 are not 
sharp. This shows the advantage of combining an £^ potential and a total variation potential. 



4.3 Experiment 3 

We revisit via the variational formulation (1.1) a pulse shape design problem investigated in [23] in 
a more restrictive setting (see also [40] for the original two-constraint formulation). This problem 
illustrates further ramifications of the proposed algorithm. 

The problem is to design a pulse shape for digital communications. The signal space is the 
standard Euclidean space Ti. = M^, where A^ = 1024 is the number of samples of the discrete pulse 
(the underlying sampling rate is 2560 Hz). Five constraints arise from engineering specifications. 
We denote by x = {S,k)o<k<N-i a signal in TC and by x = {xk)o<k<N-i its discrete Fourier 
transform. 

• The Fourier transform of the pulse should vanish at the zero frequency and at integer mul- 
tiples of 50 Hz. This constraint is associated with the set 

Ci = {x en\xln,= 0}, (4.35) 

where Di is the set of discrete frequencies at which x should vanish. 

• The modulus of the Fourier transform of the pulse should no exceed a prescribed bound 
p > beyond 300 Hz. This constraint is associated with the set 

C2 = {xen\{ykeB2)\xk\<p}, (4.36) 

where B2 represents frequencies beyond 300 Hz. 

• The energy of the pulse should not exceed a prescribed bound /z^ > in order not to interfere 
with other systems. The associated set is 

C3 = {xen\ \\x\\ <fi}. (4.37) 

• The pulse should be symmetric about its mid-point, where its value should be equal to 1. 
This corresponds to the set 

Ci = {xGn\ iN/2 = 1 and (VA; G {0, . . . , N/2}) ^fe = iN-i-k]- (4.38) 
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• The duration of the pulse should be 50 ms and it should have periodic zero crossings every 
3.125 ms. This leads to the set 

C5 = {x€ W I xl§ = 0}, (4.39) 

where S is the set of time indices in the zero areas. 

In this problem, Ci, C2, and C3 are hard constraints that must be satisfied, whereas the other 
constraints are soft ones that are incorporated via powers of distance potentials. This leads to the 
variational formulation 

minimize (£^ ix) + dP^ (x), (4.40) 

where p4 and ps are in [l,+oo[. The design problem is thus cast in the general form of (1.1), 
with 7n = 5, fi = iCi for i € {1,2,3}, and /» = dF^, for i G {4,5}. Since C3 is bounded, condi- 
tion (i) in Theorem 3.4 holds. In addition, it follows from Proposition 3.6(vi) that condition (ii) 
in Theorem 3.4 is satisfied. Indeed, 

5 
G Ci n {x G H I (Vfc G D2) IXfcl <p]r\{xen\ \\x\\ < /x} = p|ridom/i. (4.41) 

Let us emphasize that our approach is applicable to any value of (^4,^5) G [1, -l-cxop. The proximity 
operators of /4 and f^ are supplied by Proposition 2.8, whereas the other proximity operators 
are the projectors onto (Ci)i<i<3, which are straightforward [23]. A solution to (4.40) when 
P4 = P5 = 2, /9 = 10^'^'^, and ^ = 2 is shown in figure 9 and its Fourier transform is shown 
in figure 10. As is apparent in figure 9, the constraints corresponding to C4 and C5 are not 
satisfied. Forcing C4 fl C5 as a hard constraint would therefore result in an infeasible problem. 
Finally, figure 10 shows that C2 induces a 30 dB attenuation in the stop-band (beyond 300 Hz), 
in agreement with the value chosen for p. 



5 Concluding remarks 

We have proposed a proximal method for solving inverse problems that can be decomposed into 
the minimization of a sum of lower semicontinuous convex potentials. The algorithms currently 
in use in inverse problems are restricted to at most two nonsmooth potentials, which excludes 
many important scenarios and offers limited flexibility in terms of numerical implementation. By 
contrast, the algorithm proposed in the paper can handle an arbitrary number of nonsmooth 
potentials. It involves each potential by means of its own proximity operator, and activates these 
operators in parallel at each iteration. The versatility of the method is demonstrated through 
applications in signal and image recovery that illustrate various decomposition schemes, including 
one in which total variation is mixed up with other nonsmooth potentials. 
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Figure 9: Experiment 3. Pulse (amplitude versus time in ms) synthesized using 100 iterations of 
Algorithm 3.1 with 7 = 1/5. 
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Figure 10: Experiment 3. Fourier transform (amplitude in dB versus frequency in Hz) of the pulse 
of figure 9. 
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